Стабилизация перевёрнутого маятника с подвижной осью подвеса

Расчётно-графическая работа для лабораторного практикума. Вспомогательные и дополнительные материалы

Автор: В. А. Костин

2025 год $\newcommand{\ctg}{\mathop{\mathrm{ctg}}\nolimits}$ $\newcommand{\tg}{\mathop{\mathrm{tg}}\nolimits}$ $\newcommand{\arctg}{\mathop{\mathrm{arctg}}\nolimits}$ $\newcommand{\degree}{^{\circ}}$ $\renewcommand{\Re}{\mathop{\mathrm{Re}}\nolimits}$

Приложение 1. Введение в конструкции языка Python¶

Функция print и вывод ячеек¶

При настройках блокнота по умолчанию текстовое представление выходного результата последней команды (кроме инструкций присвоения) в ячейке печатается в вывод ячейки (cell output). Настройки может быть изменены установкой значения InteractiveShell.ast_node_interactivity, InteractiveShell определён в модуле IPython.core.interactiveshell. Возможные значения InteractiveShell.ast_node_interactivity: 'all', 'last', 'last_expr' (по умолчанию), 'none', 'last_expr_or_assign'.

In [5]:
y = 5 + 3
2 + 8
5 + 2
Out[5]:
7

Точка с запятой подавляет печать выходного результата последней инструкции.

In [7]:
6;

В Python 3 функция print в Python 3 по умолчанию выводит в стандартный вывод объекты, являющиеся её начальными неименованными аргументами. Стандартный вывод интерпретатора для инструкций в некоторой ячейке с кодом (code cell) в блокноте Jupyter печатается непосредственно после этой ячейки.

In [9]:
print(1, 'text', 3 + 4)

1 text 7

Форматный вывод может осуществляться с помощью метода format переменных и литералов строкового типа.

In [11]:
print('Одна треть равна {0:.3g}, а одна седьмая — {1:.6e}'.format(1 / 3,
                                                                  1 / 7))

Одна треть равна 0.333, а одна седьмая — 1.428571e-01

Арифметические операции¶

Базовые арифметические операции, в целом обозначаются так же, как и в большинстве языков. Круглые скобки могут быть использованы для группировки действий. В частности, сложение выполняется с помощью оператора +.

In [14]:
5 + 2
Out[14]:
7

Вычитание выполняется с помощью оператора -.

In [16]:
5 - 2
Out[16]:
3

Умножение выполняется с помощью оператора *.

In [18]:
5 * 2
Out[18]:
10

Вещественное деление выполняется с помощью оператора /. Результат вещественного деления — всегда число с плавающей точкой.

In [20]:
5 / 2
Out[20]:
2.5

Целочисленное деление выполняется с помощью оператора //. Результат — целое число, не превышающее вещественное частное.

In [22]:
5 // 2
Out[22]:
2

Остаток от деления может быть получен с помощью оператора %.

In [24]:
5 % 2
Out[24]:
1

Оператор ** позволяет выполнять возведение в степень.

In [26]:
5 ** 2
Out[26]:
25

Для всех операторов, перечисленных выше, имеются варианты с присвоением: +=, -=, *=. /=, //=. %=, **=.

In [28]:
y = 5
y += 2
y
Out[28]:
7

Имеется встроенная поддержка комплексной арифметики.

In [30]:
(2 + 5j) * (2 - 5j)
Out[30]:
(29+0j)

Типизация¶

Python — язык с динамической типизацией и развитой системой преобразований типов по умолчанию. Тип переменной может быть определён с помощью встроенной функции type.

In [33]:
x = 5 / 2
type(x)
Out[33]:
float
In [34]:
x = 5 // 2
type(x)
Out[34]:
int
In [35]:
x = 'text'
type(x)
Out[35]:
str

Python — строго типизованный язык.

In [37]:
x = 'text'
y = 'something'
try:
    y = x + 5
except TypeError:
    print('Нельзя прибавить число к строке!')
y

Нельзя прибавить число к строке!
Out[37]:
'something'

Операторы контроля потока выполнения¶

Python предлагает традиционный набор операторов для контроля потока выполнения. Среди них условный оператор if. Альтернативный блок указывается после оператора else. В случае нескольких условий и альтернатив может быть использован оператор elif. Имеется традиционный набор операторов сравнения ==, <, >, <=, >=, !=.

In [40]:
x = 2
if x == 3:
    print('x = 3')
elif x == 2:
    print('x = 2')
elif x == 1:
    print('x = 1')
else:
    print('x is neither 1, 2, or 3')

x = 2

Циклы могут быть организованы с помощью операторов while и for. Внутри цикла могут быть использованы операторы break (для прерывания цикла) и continue (для перехода к следующей итерации).

In [42]:
x = 10
y = x
factorial = 1
while x > 0:
    factorial *= x
    x -= 1
print('Факториал {0} равен {1}.'.format(y, factorial))

Факториал 10 равен 3628800.

Оператор цикла for используется в паре с оператором in, после которого указывается итерируемый объект (iterable). В качестве такого объекта могут выступать списки, кортежи и словари, а также некоторые другие объекты. Для итерации по арифметической прогрессии может быть использован итератор, конструируемый функцией range. Если range используется с одним целочисленным аргументом, например, range(stop), то итерация проводится от нуля последовательно (через единицу) до stop - 1. Если аргументов 2, range(start, stop), то итерация проводится от start до stop - 1. Третий аргумент может быть использован, для того чтобы указать шаг арифметической прогрессии.

In [44]:
a = 5
b = 15
sum = 0
for i in range(a, b + 1):
    sum += i
print('Сумма {0} + {1} + ... + {2} равна {3}.'.format(a, a + 1, b, sum))

Сумма 5 + 6 + ... + 15 равна 110.

Функции могут быть определены с помощью оператора def. Возвращаемый результат может быть указан с помощью оператора return. Функция может возращать несколько значений (то есть возращать кортеж значений).

In [46]:
def f1(x):
    return x + 1, 2 * x

f1(3)
Out[46]:
(4, 6)

Короткие функции могут быть определены с помощью оператора lambda.

In [48]:
f2 = lambda x: (x + 1, 2 * x)

f2(3)
Out[48]:
(4, 6)

Кортежи и списки¶

Python обладает несколькими встроенными типа составных объектов. В частности, такими типами являются списки и кортежи (неизменяемые списки). Списки конструируются из элементов с помощью квадратных скобок [], кортежи — с помощью круглых ().

In [51]:
x = [1, 2, 5]
type(x)
Out[51]:
list
In [52]:
x = (1, 2, 5)
type(x)
Out[52]:
tuple

Удобным методом создания списков и кортежей является списковое включение (list comprehension), осуществляемое с помощью операторов квадратных скобок и операторов for и in.

In [54]:
x = [1, 2, 5]
y = [i + 1 for i in x]
y
Out[54]:
[2, 3, 6]

Обращение к элементу осуществляется посредством оператора квадратных скобок [].

In [56]:
x = [1, 2, 5]
y = (1, 2, 5)
x[2], y[2]
Out[56]:
(5, 5)

Списки и кортежи в Python поддерживают срезы (slicing), то есть обращение к подсписку или подкортежу, заданному арифметической прогрессией индексов.

In [58]:
x = [i for i in range(10)]
print('Исходный список: {0}'.format(x))
print('Его слайс (срез) [3:9:2]: {0}'.format(x[3:9:2]))

Исходный список: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
Его слайс (срез) [3:9:2]: [3, 5, 7]
In [59]:
x = [i for i in range(10)]
print('Исходный список: {0}'.format(x))
print('Список в противоположном порядке: {0}'.format(x[::-1]))

Исходный список: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
Список в противоположном порядке: [9, 8, 7, 6, 5, 4, 3, 2, 1, 0]

Изменяемые и неизменяемые объекты. Сущность присваиваиния в Python¶

Любой объект в Python может быть изменяемым (mutable) или неизменяемым (immutable). Числа и строки — примеры неизменяемых объектов. Неизменяемыми являются объекты числовых типов, булевы значения, строки, неизменяемые последовательности байтов (тип bytes), кортежи (тип tuple), неизменяемые (замороженные) множества (тип frozenset). Остальные объекты, например, списки, являются изменяемыми объектами. Присваивание переменной изменяемого или неизменяемого объекта носит различный характер.

In [62]:
# Числа — неизменяемые (immutable) объекты
x = 4  # Имя «x» связывается с неизменяемым объектом «4»
y = x  # Имя «y» связывается с неизменяемым объектом «4»
y += 1 # Имя «y» связывается с неизменяемым объектом «5», значение x не
       # меняется
x
Out[62]:
4
In [63]:
# Списки — изменяемые (mutable) объекты
# объектами
x = [4]   # Выделяется память, в которую записывается список «[4]», имя «x»
          # связывается с адресом списка в памяти
y = x     # Имя «y» связывается с тем же адресом, копии объекта не создаётся
y[0] += 1 # y[0] и x[0] — синонимы
x
Out[63]:
[5]

Приложение 2. Примеры использования библиотеки Numpy¶

Подключение библиотеки¶

При работе с Python в неинтерактивном режиме для библиотеки Numpy рекомедуется использовать стандартную инструкцию импорта import numpy или для сокращения import numpy as np и обращаться к фукнциям и классам библиотеки через модуль, например, numpy.arange или np.arange. При интерактивной работе, например, в Jupyter, для подключения объектов Numpy и Matplotlib удобно воспользоваться модулем pylab и импортировать все объекты в глобальное пространство имён. При этом функции и классы оказываются доступны без указания библиотеки, то есть используется arange вместо numpy.arange.

In [67]:
from pylab import *

Многомерные массивы¶

Библиотека вводит новый тип составного объекта — numpy.ndarray. Этот тип позволяет хранить и проводить операции с многомерными массивами однородных (то есть однотипных) данных. Функция array позволяет создать такой массив из других составных объектов, например, из простых и вложенных списков или кортежей.

In [70]:
x = array([[0, 1, 2], [3, 4, 5]])
x
Out[70]:
array([[0, 1, 2],
       [3, 4, 5]])
In [71]:
type(x)
Out[71]:
numpy.ndarray

Тип элемента массива содержится в поле dtype.

In [73]:
x.dtype
Out[73]:
dtype('int32')

Поле ndim содержит число измерений массива.

In [75]:
x.ndim
Out[75]:
2

Поле size содержит число элементов массива.

In [77]:
x.size
Out[77]:
6

Поле shape содержит кортеж с размерами массива вдоль различных измерений

In [79]:
x.shape
Out[79]:
(2, 3)

Доступ к отдельному элементу может быть осуществлён посредством оператора квадратных скобок [].

In [81]:
x[1, 2], x[1][2]
Out[81]:
(5, 5)

Создание специальных массивов¶

Массив содержащий арифметическую прогрессию может быть создан с помощью функции arange. Смысл первых трёх аргументов функции аналогичен смыслу аргументов функции range, но в качестве аргументов могут выступать не только целые числа, но и числа с плавающей точкой.

In [84]:
arange(1., 4., 0.3)
Out[84]:
array([1. , 1.3, 1.6, 1.9, 2.2, 2.5, 2.8, 3.1, 3.4, 3.7])

Функция linspace также как функция arange выдаёт массив, содержащий арифметическую прогрессию, но принимает в качестве аргументов начальное значение, конечное значение и размер выходного массива.

In [86]:
linspace(2, 4, 5)
Out[86]:
array([2. , 2.5, 3. , 3.5, 4. ])

Функция meshgrid конструирует координатные матрицы из координатных векторов.

In [88]:
x = arange(1, 4)
y = arange(1, 4)
X, Y = meshgrid(x, y)
print('Матрица X:\n{0}'.format(X))
print('Матрица Y:\n{0}'.format(Y))

Матрица X:
[[1 2 3]
 [1 2 3]
 [1 2 3]]
Матрица Y:
[[1 1 1]
 [2 2 2]
 [3 3 3]]

Функция zeros создаёт массив заданных формы и типа, заполненный нулями.

In [90]:
zeros((2, 3), dtype=double)
Out[90]:
array([[0., 0., 0.],
       [0., 0., 0.]])

Функция full создаёт массив с одинаковыми значениями всех элементов.

In [92]:
full((2, 3), 5.)
Out[92]:
array([[5., 5., 5.],
       [5., 5., 5.]])

Функция empty создаёт пустой массив (то есть при его создания память не перезаписывается какими-либо заданными значениями).

In [94]:
empty((2, 3), dtype=double)
Out[94]:
array([[5., 5., 5.],
       [5., 5., 5.]])

Функции zeros_like, full_like и empty_like действуют как функции zeros, full и empty, но принимают в качестве аргумента вместо формы другой массив у которого копирует форму и тип (если тот не указывается явно среди аргументов).

In [96]:
x = zeros((2, 3), dtype=int)
empty_like(x)
Out[96]:
array([[  74515089,  913440851, 1912602624],
       [        54,      70002,    3568128]])

Так как массивы Numpy являются изменяемыми объектами, присваивание не создаёт копии массива. Для создания копии массива используется функция copy.

In [98]:
x = arange(10)
y = x
y[5] = -100
x
Out[98]:
array([   0,    1,    2,    3,    4, -100,    6,    7,    8,    9])
In [99]:
x = arange(10)
y = copy(x)
y[5] = -100
x
Out[99]:
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])

Поэлементные операции¶

Для массивов Numpy переопределены арифметические операции и элементарные математические функции. В большистве случаев операции и функции применяются поэлементно. При этом для бинарных операций оба операнда должны иметь одинаковую форму либо один из операндов должен быть скаляром. При выполнении операций тип элемента массива может меняться (например, при делении целочисленных массивов чисел или вычислении трансцендентых функций от таких массивов).

In [102]:
x = array([[0, 1, 2], [3, 4, 5]])
x
Out[102]:
array([[0, 1, 2],
       [3, 4, 5]])
In [103]:
y = x + 1
y
Out[103]:
array([[1, 2, 3],
       [4, 5, 6]])
In [104]:
x * y
Out[104]:
array([[ 0,  2,  6],
       [12, 20, 30]])
In [105]:
z = sin(x)
z
Out[105]:
array([[ 0.        ,  0.84147098,  0.90929743],
       [ 0.14112001, -0.7568025 , -0.95892427]])
In [106]:
x.dtype, z.dtype
Out[106]:
(dtype('int32'), dtype('float64'))

Срезы¶

Массивы Numpy поддерживают срезы (слайсинг) для операций с подмассивами.

In [109]:
x = arange(10)
x
Out[109]:
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])

Для положительных c срез a:b:c выбирает каждый c-й элемент с номерами ниже b, начиная с элемента номером a, при этом значения a и b берутся по модулю размера массива в соответствующем измерении.

In [111]:
x[3:9:2]
Out[111]:
array([3, 5, 7])

Для отрицательных c срез a:b:c выбирает каждый c-й элемент в обратном порядке с номерами выше b, начиная с элемента номером a.

In [113]:
x[9:3:-2]
Out[113]:
array([9, 7, 5])

Каждое из числе a, b, c в срезе a:b:c может быть опущено. Значение a и c по умолчанию равны 0 и 1 соответственно, значение b по умолчанию равно размеру массива по соответствующему измерению. Двоеточие в конце среза также может опускаться.

In [115]:
x[::-1]
Out[115]:
array([9, 8, 7, 6, 5, 4, 3, 2, 1, 0])

Одиночное двоеточие обозначает все элементы в каком-то измерении. Его удобно использовать, чтобы выбирать проекции многомерного массива, например, выбирать некоторые столбцы в двумерном массиве (матрице).

In [117]:
y = zeros((5, 7))
y[::2, ::2] = 1
y
Out[117]:
array([[1., 0., 1., 0., 1., 0., 1.],
       [0., 0., 0., 0., 0., 0., 0.],
       [1., 0., 1., 0., 1., 0., 1.],
       [0., 0., 0., 0., 0., 0., 0.],
       [1., 0., 1., 0., 1., 0., 1.]])
In [118]:
y[:, 2]
Out[118]:
array([1., 0., 1., 0., 1.])

При аналогичном выборе строк в матрице запятую с двоеточием в конце можно опустить.

In [120]:
y[2]
Out[120]:
array([1., 0., 1., 0., 1., 0., 1.])

Приложение 3. Примеры использования библиотеки Matplotlib¶

Подключение и настройка библиотеки¶

Библиотека Matplotlib предлагает два подхода к построению рисунков императивный (схожий с Matlab) и объектно-ориентированный. Для простой графики, по всей видимости, более удобен императивный подход, реализуемый с помощью функций модуля matplotlib.pyplot. При работе с Python в неинтерактивном режиме для его подключения рекомедуется использовать стандартную инструкцию импорта import matplotlib.pyplot as plt и обращаться к фукнциям и классам библиотеки через модуль, например, plt.plot. При интерактивной работе, например, в Jupyter, для подключения объектов Numpy и Matplotlib удобно воспользоваться модулем pylab и импортировать все объекты в глобальное пространство имён. При этом функции и классы оказываются доступны без указания библиотеки, то есть используется plot вместо plt.plot.

In [124]:
from pylab import *

Так называемая «магическая» команда %matplotlib позволяет выбрать тип вывода (бэкенд) рисунков. Параметр inline указывает на то, что рисунки Matplotlib должны быть встроены как растровые изображения, для интерактивных рисунков на JavaScript можно использовать параметр notebook вместо inline. Команда %matplotlib --list выводит список возможных типов вывода (бэкендов).

In [126]:
%matplotlib inline

Переменная rcParams содержит различные настройки рисунков. Доступ к настройкам осуществляется по ключам. Например, ключи figure.figsize и figure.dpi соответствуют размеру рисунков по умолчанию в дюймах и разрешению рисунков в точках на дюйм.

In [128]:
rcParams['figure.figsize'] = (9.6, 7.2)
rcParams['figure.dpi'] = 100

Построение графика функции¶

Один из наиболее простых способов построить график функции или параметрически заданную кривую в Matplotlib — использовать функцию plot. Функция принимает в качестве аргумента массивы абсцисс и ординат точек кривой. Рисунок можно снабдить заголовком с помощью функции title, а подписать оси можно с помощью функций xlabel и ylabel. Добавить подпись в произвольном месте можно с помощью команды text. Подписи могут содержать формулы в формате, схожем с LaTeX.

In [131]:
x = arange(0., 10., 0.1)
y = sin(x)
plot(x, y)
xlabel(r'$x$') # Подписи по осям могут содержать формулы
ylabel(r'$\sin x$') # r перед строкой позволяет не экранировать специальные 
                    # символы, в частности, обратную косую черту \
text(3., 0.5, r'Кривая 1')
title('График синуса'); # Точка с запятой подавляет вывод последней команды
No description has been provided for this image

Функция plot также может принимать аргументы, контролирующие параметры отображения кривой.

In [133]:
plot(x, y, linewidth=4.0)
xlabel(r'$x$')
ylabel(r'$\sin x$')
title('График синуса');
No description has been provided for this image

Третьим аргументов функция plot может принимать стилевую строку, указывающую цвет кривой, её тип и вид маркеров.

In [135]:
plot(x, y, 'ro')
xlabel(r'$x$')
ylabel(r'$\sin x$')
title('График синуса');
No description has been provided for this image

Функция plot может несколько наборов абсцисс, ординат и стилевых строк для построения нескольких графиков в одних осях.

In [137]:
plot(x, y, x[::9], y[::9], 'ro', x, cos(x))
xlabel(r'$x$')
ylabel(r'$\sin x$ и $\cos x$')
title('Графики синуса и косинуса');
No description has been provided for this image

Несколько графиков можно построить, вызвав функцию plot несколько раз.

In [139]:
plot(x, y)
plot(x[::9], y[::9], 'ro')
plot(x, cos(x))
xlabel(r'$x$')
ylabel(r'$\sin x$ и $\cos x$')
title('Графики синуса и косинуса');
No description has been provided for this image

Функция fill аналогична функции plot, но вместо графика или точек рисует внутреннюю область графика, заполненную цветом. В примере ниже инструкция axis('equal') устанавливает одинаковые масштабы по горизонтальной и вертикальной осям.

In [141]:
t = linspace(0., 2 * pi, 201)
fill(cos(t), sin(t), color='pink')
fill([2, 2, 4, 4], [-1, 1, 1, -1], color='purple')
axis('equal')
xlabel(r'$x$')
xlabel(r'$y$')
title('Круг и квадрат');
No description has been provided for this image

Построение нескольких координатных плоскостей (осевых коробок) на одном рисунке¶

Несколько координатных плоскостей (осей, осевых коробок, axes) на одном рисунке можно получить с помощью функции subplot. Команда subplot(nrows, ncols, index) создаёт виртуальную сетку из координатных плоскостей с nrows строками и ncols столбцами и возвращает координатную плоскость с индексом index из этой сетки. Общий заголовок для всех координатных плоскостей можно получить с помощью функции suptitle. Функция tight_layout может использовать для корректировки размеров осей тем, чтобы подписи по осям не перекрывали других подрисунков и не выходили за пределы всего рисунка. Аргумент rect указывает на положение граничной рамки, содержащей все подрисунки, но не главный заголовок.

In [144]:
x = arange(0., 10., 0.1)
subplot(2, 1, 1) # Первая панель (первый «подрисунок»)
plot(x, sin(x))
xlabel(r'$x$')
ylabel(r'$\sin x$')
subplot(2, 1, 2) # Вторая панель (второй «подрисунок»)
plot(x, cos(x))
xlabel(r'$x$')
ylabel(r'$\cos x$')
suptitle(r'Графики синуса и косинуса')
tight_layout(rect=[0, 0, 1, 0.95]);
No description has been provided for this image

Построение фазового портрета¶

Функцию streamplot удобно использовать для построения фазовых портретов, силовых линий векторных полей или линий тока на плоскости. В простейшем варианте функция принимает четыре аргумента: первые два задают координатные матрицы, два следующих задают две компоненты векторного поля. Одним из возможных дополнительных параметров является density, определяющий характерную плотность отображаемых линий (то есть обратное расстояние между линиями). Далее приводится пример построения фазовой портрета с использованием функции streamplot для системы $\ddot\varphi + 0.3\dot\varphi + \sin\varphi = 0$.

In [147]:
# Создание сеток переменных 201 на 201 элемента
phi, dphi = meshgrid(linspace(-pi, pi, 201), linspace(-3., 3., 201))

# Вычисление вертикальной компоненты векторного поля (горизонтальная
# компонента равна dphi и её не нужно специально вычислять)
ddphi = -sin(phi) - 0.3 * dphi

# Построение фазовой плоскости: первые два аргумента задают сетку значений
# абсцисс и ординат, следующие два аргумента задают горизонтальную и
# вертикальную компоненты векторного поля
streamplot(phi, dphi, dphi, ddphi, density=1.5, color='k')

xlabel(r'$\varphi$')
ylabel(r'$\dot\varphi$')
title(r'Фазовый портрет системы $\ddot\varphi + 0.3\dot\varphi + ' +
      r'\sin\,\varphi = 0$');
No description has been provided for this image

Построение контурной карты¶

Линии уровня некоторой функции двух переменных можно построить с помощью функции contour. В одном из вариантов функция принимает первыми аргументами две координатные матрицы и матрицу значений функции. Четвёртый аргумент может задать количество линий уровня или упорядоченный по возрастанию набор значений фукнции для которых строятся линии уровня. Функция clabel может быть использована для расстановки меток значений уровня на линиях. В этой функции параметр inline контролирует, стираются ли линии уровня под метками, параметр fontsize управляет размером шрифта меток. Рисунок также можно снабдить шкалой с помощью функции colorbar.

In [150]:
x = linspace(-2., 2., 200)
y = linspace(-2., 2., 200)
X, Y = meshgrid(x, y)
Z = X * exp(-X ** 2 - 1.5 * Y ** 2) + 0.2 * Y
v = linspace(-0.5, 0.5, 21)
contours = contour(X, Y, Z, v)
gca().set_aspect('equal') # Установка одинаковых масштабов по вертикальной
                          # и горизонтальной осям (здесь необязательно)
clabel(contours, inline=True, fontsize=8)
colorbar(contours, orientation='vertical', label=r'$z$')
xlabel(r'$x$')
ylabel(r'$y$')
title(r'Линии уровня $z = e^{-x^2 - 1.5y^2} + 0.2y$');
No description has been provided for this image

Функция contourf работает аналогично contour но строит не сами линии уровня, а закрашивает разными цветами области между линиями уровня. Такой рисунок также можно снабдить шкалой с помощью функции colorbar.

In [152]:
x = linspace(-2., 2., 200)
y = linspace(-2., 2., 200)
X, Y = meshgrid(x, y)
Z = X * exp(-X ** 2 - 1.5 * Y ** 2) + 0.2 * Y
v = linspace(-0.5, 0.5, 21)
contours = contourf(X, Y, Z, v)
gca().set_aspect('equal') # Установка одинаковых масштабов по вертикальной
                          # и горизонтальной осям (здесь необязательно)
colorbar(contours, orientation='vertical', label=r'$z$')
xlabel(r'$x$')
ylabel(r'$y$')
title(r'Тепловая карта $z = e^{-x^2 - 1.5y^2} + 0.2y$');
No description has been provided for this image

Построение тепловой карты¶

Упомянутой выше в разделе 5.3.5 функции contourf может использоваться для построения векторной тепловой карты, представляющей набор областей равных цветов. Для построения растровой тепловой карты можно применять функцию imshow. В этой функции параметр origin указывает, верхнему или нижнему левому углу панели соответствует первый (с индексами 0, 0) элемент изображаемого массива. Обычно этот параметр должен быть установлен в 'lower' при построении тепловых карт функций, заданных матрицами в традиционных декартовых координатах. Параметр extent указывает координаты левого нижнего и правого верхнего углов тепловой карты.

Приложение 4. Справка по функциям модуля control_theory¶

Plant¶

class Plant()

Класс Plant — является базовым интерфейсным классом для представления объекта управления, задаваемого с помощью системы обыкновенных дифференциальных уравнений $$\dot{\mathbf x} = \mathbf f(\mathbf x, \mathbf u, t),$$ где $\mathbf x$ — векторная переменная, определяющеая состояние объекта (элемент пространства состояний), $\mathbf u$ — векторное управление, функция $\mathbf f$ полностью определяет модель объекта управления.

Экземпляр базового класса Plant представляет пустой (тривиальный) объект управления без переменных, параметров и управления. Экземпляры классов-наследников представляют специфичные нетривиальные объекты управления.

Plant.PLANT_DIM, Plant.CONTROL_DIM, Plant.OTYPE¶

Класс и его наследники содержит следующие поля класса (то есть поля, общие для всех экземпляров класса).

Название Тип Описание
PLANT_DIM int Размерность состояния объекта (то есть вектора $\mathbf x$).
CONTROL_DIM int Размерность управления (то есть вектора $\mathbf u$).
OTYPE str Название (строковый идентификатор) типа объекта управления.

Plant.__call__¶

f = Plant.__call__(x, u=zeros(CONTROL_DIM), t=0.)

Экзмепляр класса является функтором и при вызове вычисляет значение временной производной $\dot {\mathbf x}$, иными словами вычисляет функцию $\mathbf f$.

Принимаемые аргументы

Название Тип Описание Значение по умолчанию
x array_like, shape=(PLANT_DIM,) Вектор состояния объекта $\mathbf x$.
u array_like, shape=(CONTROL_DIM,) Вектор управления $\mathbf u$ (по умолчанию — нулевой). zeros(CONTROL_DIM)
t number Значение времени $t$. 0.

Возвращаемые значения

Название Тип Описание
f ndarray, shape=(PLANT_DIM,), dtype='double' Вектор временной производной $\mathbf f(\mathbf x, \mathbf u, t)$ состояния объекта.

Plant.__str__¶

plant_str = Plant.__str__()

Для экзмепляра класса определена функция преобразования в строку __str__.

Возвращаемые значения

Название Тип Описание
plant_str str Строковое представление экземпляра класса.

PendulumTorque¶

Класс PendulumTorque — потомок класса Plant для представления маятника, описываемого уравнением $\ddot\varphi + \nu\dot\varphi - \Omega^2 \sin\varphi = -u$, где $(\varphi, \dot\varphi)$ — вектор состояния и $u$ — одномерное управление.

PendulumTorque.omega, PendulumTorque.nu¶

Экземпляр класса содержит следующие поля, соответствующие параметрам системы. При инициализации объекта связываются со случайными действительными числами.

Название Тип Описание
omega float Параметр $\Omega$.
nu float Параметр $\nu$.

PendulumHorizontalForce¶

Класс PendulumHorizontalForce — потомок класса Plant для представления маятника, описываемого уравнением $\ddot\varphi + \nu\dot\varphi - \Omega^2 \sin\varphi = -u\cos\varphi$, где $(\varphi, \dot\varphi)$ — вектор состояния и $u$ — одномерное управление.

PendulumHorizontalForce.omega, PendulumHorizontalForce.nu¶

Экземпляр класса содержит следующие поля, соответствующие параметрам системы $(4)$. При инициализации объекта связываются со случайными действительными числами.

Название Тип Описание
omega float Параметр $\Omega$.
nu float Параметр $\nu$.

Control¶

class Control()

Класс Control — является базовым интерфейсным классом для представления стратегии управления, задаваемой с помощью системы уравнений $$\mathbf u = \mathbf g(\mathbf x, \mathbf v, t),$$ $$\dot{\mathbf v} = \mathbf h(\mathbf x, \mathbf v, t),$$ где $\mathbf x$ — векторная переменная, определяющеая состояние объекта (элемент пространства состояний), $\mathbf u$ — векторное управление, $\mathbf v$ — внутренние («скрытые») переменные, характеризующие состояние регулятора, функция $\mathbf h$ определяет внутреннюю динамику скрытых переменных регулятора, функция $\mathbf g$ связывает управление со состояниями объекта управления и регулятора.

Экземпляр базового класса Control представляет пустой регулятор без переменных. Экземпляры классов-наследников представляют специфичные нетривиальные регуляторы.

Control.PLANT_DIM, Control.CONTROL_DIM, Control.HIDDEN_DIM, Control.CTYPE¶

Класс и его наследники содержит следующие поля класса (то есть поля, общие для всех экзмепляров класса).

Название Тип Описание
PLANT_DIM int Размерность состояния объекта (то есть вектора $\mathbf x$).
CONTROL_DIM int Размерность управления (то есть вектора $\mathbf u$).
HIDDEN_DIM int Размерность состояния регулятора (то есть вектора $\mathbf v$).
CTYPE str Название (строковый идентификатор) типа регулятора.

Control.__call__¶

g = Control.__call__(x, v, t=0.)

Экзмепляр класса является функтором и при вызове вычисляет значение временной производной $\dot {\mathbf v}$, иными словами вычисляет функцию $\mathbf h$.

Принимаемые аргументы

Название Тип Описание Значение по умолчанию
x array_like, shape=(PLANT_DIM,) Вектор состояния объекта $\mathbf x$.
v array_like, shape=(HIDDEN_DIM,) Вектор состояния регулятора $\mathbf v$.
t number Значение времени $t$. 0.

Возвращаемые значения

Название Тип Описание
g ndarray, shape=(HIDDEN_DIM,), dtype='double' Вектор временной производной $\mathbf h(\mathbf x, \mathbf v, t)$ состояния регулятора.

Control.control¶

u = Control.control(x, v, t=0.)

Метод control вычисляет управление $\mathbf u$ по состояниям объекта и регулятора, иными словами вычисляет функцию $\mathbf g$.

Принимаемые аргументы

Название Тип Описание Значение по умолчанию
x array_like, shape=(PLANT_DIM,) Вектор состояния объекта $\mathbf x$.
v array_like, shape=(HIDDEN_DIM,) Вектор состояния регулятора $\mathbf v$.
t number Значение времени $t$. 0.

Возвращаемые значения

Название Тип Описание
u ndarray, shape=(CONTROL_DIM,), dtype='double' Вектор управления $\mathbf u$.

Control.__str__¶

control_str = Control.__str__()

Для экзмепляра класса определена функция преобразования в строку __str__.

Возвращаемые значения

Название Тип Описание
control_str str Строковое представление экземпляра класса.

ZeroControl¶

class ZeroPendulumControl(plant_dim, control_dim)

Класс ZeroPendulumControl — потомок класса Control для представления нулевого управления для маятником на тележке Pendulum.

Принимаемые аргументы

Название Тип Описание Значение по умолчанию
plant_dim int Размерность состояния объекта управления. 1
control_dim int Размерность управления. 1

StateFeedbackControl¶

class StateFeedbackControl(plant_dim, control_dim, fun)

Класс StateFeedbackControl — потомок класса ZeroControl для представления управления по состоянию объекта.

Принимаемые аргументы

Название Тип Описание Значение по умолчанию
plant_dim int Размерность состояния объекта управления. 1
control_dim int Размерность управления. 1
fun u = callable(x); x: arraylike, shape=(plant_dim,); u: arraylike, (control_dim,) или None Функция, задающая стратегию управления: принимает текущее состояние объекта x и возращает управление u. Если передаётся None, то функция считается возвращающей нулевой вектор (поведение класса не отличается от ZeroControl). None

StateFeedbackControl.fun¶

Экземпляр класса содержит поле с функцией, задающей стратегию управления.

Название Тип Описание
fun u = callable(x); x: arraylike, shape=(PLANT_DIM,); u: arraylike, (CONTROL_DIM,) Функция, задающая стратегию управления: принимает текущее состояние объекта x и возращает управление u.

StateTimeFeedbackControl¶

class StateTimeFeedbackControl(plant_dim, control_dim, fun)

Класс StateFeedbackControl — потомок класса ZeroControl для представления комбинированного управления по состоянию объекта и моменту времени (программного управления).

Принимаемые аргументы

Название Тип Описание Значение по умолчанию
plant_dim int Размерность состояния объекта управления. 1
control_dim int Размерность управления. 1
fun u = callable(x, t); x: arraylike, shape=(plant_dim,); t: number; u: arraylike, (control_dim,) или None Функция, задающая стратегию управления: принимает текущее состояние объекта x и время t и возращает управление u. Если передаётся None, то функция считается возвращающей нулевой вектор (поведение класса не отличается от ZeroControl). None

StateTimeFeedbackControl.fun¶

Экземпляр класса содержит поле с функцией, задающей стратегию управления.

Название Тип Описание
fun u = callable(x, t); x: arraylike, shape=(PLANT_DIM,); t: number; u: arraylike, (CONTROL_DIM,) Функция, задающая стратегию управления: принимает текущее состояние объекта x и время t и возращает управление u.

LinearStateControl¶

class LinearStateControl(plant_dim, control_dim, k=None)

Класс LinearStateControl — потомок класса ZeroControl для представления линейного управления по состоянию.

Принимаемые аргументы

Название Тип Описание Значение по умолчанию
plant_dim int Размерность состояния объекта управления. 1
control_dim int Размерность управления. 1
k arraylike, shape=(control_dim, plant_dim) или shape=(plant_dim,), если control_dim равно 1, или None Матрица или вектор $\mathbf k$ коэффициентов линейного закона управления. Задание None эквивалентно заданию нулевой матрицы с shape=(control_dim, plant_dim). None

LinearStateControl.k¶

Экземпляр класса содержит поле с матрицей $\mathbf k$ коэффициентов линейного закона управления.

Название Тип Описание
k ndarray, shape=(CONTROL_DIM, PLANT_DIM), dtype='double' Матрица коэффциентов $\mathbf k$ линейного закона управления (матрица умножается на вектор слева).

rk¶

tout, y = rk(rhs, y0, t, N=1, stop_condition=None)

Функция rk реализует классический метод Рунге — Кутты четвёртого порядка для решения обыкновенных дифферециальных уравений.

Принимаемые аргументы

Название Тип Описание Значение по умолчанию
rhs dydt1 = callable(y1, t1); y1: ndarray, dtype=float or dtype=complex; t1: number; dydt1: ndarray, shape=y1.shape Функция, задающая систему обыкновенных уравнений: для массива y1 значений зависимых переменных (координат) и независимой переменной (времени) t1 возвращает массив dydt1 значений производных координат по времени.
y0 arraylike, shape=y1.shape Массив начальных значений координат.
t arraylike, ndim=1 Сортированный по возрастанию или убыванию массив значений времени, для которых вычисляются значения координат. Первый элемент массива отвечает времени для которого задаётся начальное значение. Если массив сортирован по убыванию, расчёт ведётся в обратном времени (то есть в сторону уменьшения времени).
N int Число шагов метода Рунге — Кутты, приходящееся на один интервал, образуемый соседними элементами из t. 1
stop_condition None или stop_flag = callable(y1, t1); y1: ndarray, dtype=float or dtype=complex; t1: number; stop_flag: bool Функция, задающая критерий остановки, если для текущих значений y1 и t1 выдаёт True, выполнение останавливается и выдаются предшествующие рассчитанные значения времени и переменного вектора. Если задано None выполнение останавливается, когда пройден весь массив времён t. None

Возращаемые значения

Название Тип Описание
tout ndarray, ndim=y1.ndim + 1, ndim=1, shape[0] <= t.shape[0] Массив, содержащий значения времён, для которых были рассчитаны значения координат. Если функция stop_condition не задана или не выдавала True в ходе расчёта, выводится копия t. Иначе выводится обрезанный вариант массива t, при этом последний элемент массива отвечает времени, на котором условие остановки было выполнено.
y ndarray, ndim=y1.ndim + 1, shape=(tout.shape[0],) + y1.shape Многомерный массив, содержащий рассчитанные значения координат для времён из tout.

pcrhs¶

rhs = pcrhs(p, c)

Функция pcrhs для заданных объекта и регулятора создаёт функцию (функциональный объект), которая рассчитывает правую часть системы обыкновенных дифференциальных уравнений, описывающих совместную эволюцию состояний объекта и регулятора во времени $$\dot{\mathbf X} = \mathbf F(\mathbf X, t),$$ где $\mathbf X = \begin{bmatrix} \mathbf x \\ \mathbf v\end{bmatrix}$ — объединённый вектор состояний объекта и регулятора, $\mathbf F(\mathbf X, t) = \begin{bmatrix} \mathbf f[\mathbf x, \mathbf g(\mathbf x, \mathbf v, t), t] \\ \mathbf h(\mathbf x, \mathbf v, t) \end{bmatrix}$. Эта система эквивалентна системе приведённых ранее уравнений $$\dot{\mathbf x} = \mathbf f(\mathbf x, \mathbf u, t),$$ $$\mathbf u = \mathbf g(\mathbf x, \mathbf v, t),$$ $$\dot{\mathbf v} = \mathbf h(\mathbf x, \mathbf v, t).$$

Принимаемые аргументы

Название Тип Описание Значение по умолчанию
p Plant Объект управления.
c Control, PLANT_DIM=p.PLANT_DIM, CONTROL_DIM=p.CONTROL_DIM Регулятор.

Возращаемые значения

Название Тип Описание
rhs dydt1 = callable(y1, t1); y1: ndarray, dtype=float or dtype=complex, shape=(p.PLANT_DIM + c.HIDDEN_DIM,); t1: number; dydt1: ndarray, shape=y1.shape Функция, вычисляющая зависимость $\mathbf F(\mathbf X, t)$.

control_output¶

ufun = control_output(c)

Функция control для заданного регулятора создаёт функцию (функциональный объект), которая рассчитывает управление по состояниям объекта и управления, $\mathbf u(\mathbf X, t) = \mathbf g(\mathbf x, \mathbf v, t)$.

Принимаемые аргументы

Название Тип Описание Значение по умолчанию
c Control, PLANT_DIM=p.PLANT_DIM, CONTROL_DIM=p.CONTROL_DIM Регулятор.

Возращаемые значения

Название Тип Описание
ufun u = callable(y1, t1); y1: ndarray, dtype=float or dtype=complex, shape=(c.PLANT_DIM + c.HIDDEN_DIM,); t1: number; u: ndarray, shape=(c.CONTROL_DIM,) Функция, вычисляющая зависимость $\mathbf u(\mathbf X, t)$.

integrate¶

t, y = integrate(p, c, x0, v0, dt, T, N=1, method=rk, return_control=False, stop_condition=None)

или

t, y, u = integrate(p, c, x0, v0, dt, T, N=1, method=rk, return_control=False, stop_condition=None)

Функция integrate для заданных объекта и регулятора рассчитывает эволюцию их состояний во времени посредством численного интегрирования соответствующих обыкновенных дифференциальных уравнений.

Принимаемые аргументы

Название Тип Описание Значение по умолчанию
p Plant Объект управления.
c Control, PLANT_DIM=p.PLANT_DIM, CONTROL_DIM=p.CONTROL_DIM Регулятор.
x0 array_like, shape=(p.PLANT_DIM,) Начальное состояние объекта.
v0 array_like, shape=(c.HIDDEN_DIM,) Начальное состояние регулятора.
dt number Шаг по времени в выходных массивах. Если задаётся отрицательным, то расчёт ведётся в обратном времени.
T number Максимальное (для положительного dt) или минимальное (для отрицательного dt) время, до которого рассчитывается эволюция.
N int Число шагов метода интегрирования, приходящееся на один шаг по времени в выходных массивах. 1
method callable Функция, реализующая метод интегрирования с сигнатурой, как у rk. rk
return_control Если True, то возвращает массив значений управления помимо массивов времени и состояний. False
stop_condition None or b1 = callable(y1, t1); y1: ndarray, shape=(p.PLANT_DIM + c.HIDDEN_DIM,); t1: number; b1: bool Необязательная функция координат и времени, задающая условие остановки расчёта, если функция возвращает True после шага, расчёт останавливается. None

Возращаемые значения

Название Тип Описание
t ndarray, ndim=1 Эквидистантный массив значений времени, для которых выводятся состояния объекта и регулятора. Если функция stop_condition не задана или не сработала в ходе расчёта, то форма выходного массива shape=(int(ceil(T/dt)),). Если функция stop_condition сработала, то последний элемент массива отвечает времени, на котором условие остановки было выполнено.
y ndarray, shape=(t.shape[0], p.PLANT_DIM + c.HIDDEN_DIM) Двумерный массив, содержащий рассчитанные значения состояний объекта и регулятора для времён из t. Состояния объекта содержатся в первых p.PLANT_DIM столбцах. В оставшихся c.HIDDEN_DIM столбцах содержатся состояния регулятора.
u ndarray, shape=(t.shape[0], p.CONTROL_DIM) Двумерный массив, содержащий рассчитанные значения управления для времён из t. Возвращается, если return_control задано как True.

animate_pendulum¶

anim = animate_pendulum(t, state, phi_lims=(-numpy.pi, numpy.pi), cylinder_mode=True, phase_portrait=None, resolution=(960, 540), dpi=108, spacing=1, invsec=1.0, filename=None, codec=None, progress=True)

Функция animate_pendulum строит объект-анимацию matplotlib.animation.Animation по заданным состояниям маятника на эквидистантной временной сетке. Анимация дополняется соответствующей траекторией на фазовом портрете на плоскости $(\varphi, \dot\varphi)$ или соответствующей развёртке фазового цилиндра.

Принимаемые аргументы

Название Тип Описание Значение по умолчанию
t ndarray, ndim=1 Сортированный по возрастанию эквидистантный массив значений времени, для которых задаются состояния маятника.
state ndarray, ndim=2, shape[0]=t.shape[0], shape[1]>=2 Двумерный массив, содержащий состояния маятника для времён из t. По первому измерению меняется время, по второму — номер переменной, фактически используются только первый и второй столбцы, в которых содержатся значения угла $\varphi$ и угловой скорости $\dot\varphi$ соответственно.
phi_lims arraylike, shape=(2,) Горизонтальные пределы фазового потрета. (-numpy.pi, numpy.pi)
cylinder_mode bool Если задано как True, то фазовый портрет предполагается развёрткой цилиндра и углы (абсциссы), отличающиеся на величины кратные $2\pi$ полагаются идентичными, а фазовый портрет является периодическим по горизонтали с периодом $2\pi$. True
phase_portrait None or (Plant, Control), Plant.PLANT_DIM=2 Если задано как None, то на фазовом портрете изображается только одна траектория, определяемая содержим state. Если задаётся пара двумерного объекта управления и самого регулятора, то рисуются дополнительные траектории, определяемые правыми частями замкнутой системы (скрытые переменные регулятора полагаются нулями). None
resolution arraylike, shape=(2,) Разрешение кадра. (960, 540)
dpi number Число пикселей на дюйм (влияет на толщину линий и размер символов). 108
spacing int Количество шагов временной сетки t на один кадр анимации. 1
invsec number Продолжительность единицы времени (в которой задана сетка t) в секундах. 1.0
filename str or None Имя файла, в который записывается видео анимации, например 'test.mp4'. Если filename задано как None, запись в файл не производится. None
codec str or None Кодировщик видео, например 'h264'. Если задано как None используется кодировщик по умолчанию. None
progress bool Если задано как True, при сохранение в файл показывается индикатор прогресса. True

Возращаемые значения

Название Тип Описание
anim matplotlib.animation.Animation Анимированный объект matplotlib.

axis_pi_ticks¶

axis_pi_ticks(axis, major_period=numpy.pi/2, minor_period=numpy.pi/6, denominator=2)

Функция axis_pi_ticks изменяет функции локатора и форматтера рисок вдоль оси axis, чтобы периодичность основных рисок составляла major, а дополнительных — minor. Форматтер форматирует метки рисок, полагая их кратными np.pi/denominator.

Принимаемые аргументы

Название Тип Описание Значение по умолчанию
axis matplotlib.axis.Axis Ось рисунка matplotlib.
major_period number Период основных рисок. numpy.pi/2
minor_period number Период дополнительных рисок. numpy.pi/6
denominator int Знаменатель $n$ дроби $\pi/n$, определяющей дискретизацию для форматтера. 2

Приложение 5. Вспомогательные примеры к работе¶

Инициализация демонстрационного маятника¶

In [210]:
from control_theory import *
p = PendulumTorque(1, 0.3)

Пример расчёта временной эволюции демонстрационного маятника с линейным управлением¶

In [212]:
c = LinearStateControl(2, 1, [2, 0])
t, s, u = integrate(p, c, x0=[0., 2], v0=empty(0), dt=0.01, T=30,
                    N=1, return_control=True)
Пример построения графика зависимости угла поворота маятника от времени¶
In [214]:
plot(t, s[:, 0])
xlabel(r'Время $t$')
ylabel(r'Угол $\varphi$')
title('Зависимость угла поворота от времени');
No description has been provided for this image
Пример построения временных зависимостей угла поворота управления на одном рисунке¶
In [216]:
subplot(2, 1, 1) # Первая панель (первый «подрисунок»)
plot(t, s[:, 0])
xlabel(r'Время $t$')
ylabel(r'Угол $\varphi$')
subplot(2, 1, 2) # Вторая панель (второй «подрисунок»)
plot(t, u)
xlabel(r'Время $t$')
ylabel(r'Управление $u$')
suptitle(r'Угол и управление')
tight_layout(rect=[0, 0, 1, 0.95]); # Корректировка размеров осей с тем, чтобы
                                    # подписи по осям не перекрывали других
                                    # подрисунков и не выходили за пределы
                                    # всего рисунка. Аргумент rect указывает
                                    # на положение граничной рамки, содержащей
                                    # все подрисунки, но не главный заголовок.
No description has been provided for this image
Пример построения фазовой траектории¶
In [218]:
plot(s[:, 0], s[:, 1])
gca().set_aspect('equal') # Установка одинаковых масштабов по вертикальной
                          # и горизонтальной осям
xlabel(r'$\varphi$')
ylabel(r'$\dot\varphi$')
title('Фазовая траектория');
No description has been provided for this image
Пример построения анимации движения маятника с выводом в блокнот и видео-файл¶
In [220]:
# Создание анимации и запись в файл `test.mp4`
anim = animate_pendulum(t, s, phi_lims=(-0.5 * pi, 0.5 * pi), cylinder_mode=False, spacing=3, invsec=1,
                         filename='files/test.mp4')

files/test.mp4, saving progress
In [221]:
# Импорт функции `HTML` для отображения произвольного кода HTML в выводе
# ячейки
from IPython.display import HTML
In [222]:
HTML(anim.to_html5_video()) # Вывод анимации в блокнот как встроенное видео
                             # HTML 5
Out[222]:
Your browser does not support the video tag.
In [223]:
rcParams['animation.embed_limit'] = 5 # Максимальный размер анимации в МБ
In [224]:
# Создание анимации меньщего размера
anim_small = animate_pendulum(t, s, phi_lims=(-0.5 * pi, 0.5 * pi), cylinder_mode=False, spacing=30, invsec=1)

# Вывод анимации в виде кода HTML
HTML(anim_small.to_jshtml())
Out[224]:
No description has been provided for this image

В качестве альтернативы можно выводить видео из предварительно созданного файла в ячейках Markdown. В таком случае видео будет отображаться и после перезагрузки страницы. Видео ниже вставлено с помощью команды ниже.

<video controls frameborder="0" allowfullscreen autoplay src="files/test.mp4" type=video/mp4></video>

Пример построения границы области притяжения положения равновесия $(0, 0)$ системы $\dot x = y$, $\dot y = -x - y + x^3$¶

В системе $\dot x = y$, $\dot y = -x - y + x^3$ имеется две седловые точки $(x, y) = (\pm 1, 0)$. Система симметрична относительно замены $(x, y) \rightarrow (-x, -y)$, и можно рассмотреть только седло $(x, y) = (1, 0)$ и две его устойчивые сепаратрисы. Вблизи этого седла система линеаризуется как $\dot x = y$, $\dot y = (x - 1) - y$ или $$\frac{d}{dt} \begin{bmatrix} x - 1 \\ y \end{bmatrix} = \begin{bmatrix} 0 & 1 \\ 1 & -1 \end{bmatrix} \begin{bmatrix} x - 1 \\ y \end{bmatrix}.$$ Матрица системы имеет собственные значения $\dfrac{-1 \pm \sqrt{5}}{2}$ и собственные векторы $\begin{bmatrix} 1 \\ \frac{-1 \pm \sqrt{5}}{2} \end{bmatrix}$. Устойчивым сепаратрисам отвечает отрицательное собственное число (нижний знак). Таким образом, в качестве начальных условий для нахождения сепаратрис можно взять $$ \begin{bmatrix} x(0) - 1 \\ y(0) \end{bmatrix} = \pm\delta \begin{bmatrix} 1 \\ \frac{1 - \sqrt{5}}{2} \end{bmatrix}, $$ где $\delta$ — некоторое достаточное малое (в масштабах нелинейности системы, то есть единицы) положительное число.

In [228]:
delta = 1e-3 # Малое число

# Пределы по осям на рисунке
xmin = -3.
xmax = 3.
ymin = -3.
ymax = 3.

# Массив времён для сепаратрисы, приходящей в точку (1, 0) снизу справа
t1 = arange(0., -8.5, -0.01)
# Массив времён для сепаратрисы, приходящей в точку (1, 0) сверху слева
t2 = arange(0., -9.5, -0.01)

# Функция вычисляет правую часть системы (векторное поле потока)
def rhs_example(y, t):
    return array([y[1], -y[1] - y[0] * (1. - y[0] * y[0])])

# Функция задаёт критерий остановки расчёта сепаратрис: расчёт
# останавливается, когда сепратриса покидает заданные границы осевой коробки
# (пределы по осям). В общем случае сепартриса, покинув заданные пределы,
# может затем вернуться обратно и могут потребоваться более сложные критерии
# остановки
def stop_condition_example(y, t):
    return y[0] > xmax or y[0] < xmin or y[1] > ymax or y[1] < ymin

# Интегрирование методом Рунге — Кутты вдоль сепаратрис системы, определяемой
# функцией rhs_example в обратном времени. Для маятника можно вместо функции
# rk использовать функцию integrate, при этом нужно указать отрицательные
# значения шага по времени dt и времени расчёта T.
_, y1 = rk(rhs_example, [1. + delta, -0.5 * (1. + sqrt(5.)) * delta], t1, 10,
           stop_condition_example)
_, y2 = rk(rhs_example, [1. - delta, 0.5 * (1. + sqrt(5.)) * delta], t2, 10,
           stop_condition_example)

# Объединение двух сепаратрис в одну кривую
y = vstack([y1[::-1], y2])
# Объединение двух симметричных пар сепаратрис в одну кривую
y = vstack([y, -y])

# Заполнение области, ограниченной сепаратрисами цветом
fill(y[:, 0], y[:, 1], color='pink')
# Построение сепаратрис
plot(y[:, 0], y[:, 1], color='red')
# Точки в положениях равновесия
plot([0], [0], 'bo')
plot([-1, 1], [0, 0], 'ro')

# Установление пределов по осям
axis([xmin, xmax, ymin, ymax])

# Построение фазовой плоскости
x, y = meshgrid(linspace(xmin, xmax, 201), linspace(ymin, ymax, 201))
dy = -y - x * (1. - x * x)
streamplot(x, y, y, dy, color='k', density=1.5)

xlabel(r'$x$')
ylabel(r'$y$')
title(r'Область притяжения устойчивого положения равновесия системы ' +
      r'$\dot x = y$, $\dot y = -x - y + x^3$');
No description has been provided for this image